{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 背景\n",
    "现实世界中的时间序列预测具有很大的挑战性，其原因不仅限于问题特征，例如具有多个输入变量，需要预测多个时间步长以及需要对多个物理站点执行相同类型的预测。\n",
    "\n",
    "EMC数据科学全球Hackathon数据集，或简称为“空气质量预测”数据集，描述了多个站点的天气状况，并要求对随后三天的空气质量测量值进行预测。\n",
    "\n",
    "机器学习算法可以应用于时间序列预测问题，并提供诸如处理带有复杂杂讯的多个输入变量的功能。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 问题描述\n",
    "空气质量预测数据集描述了多个地点的天气状况，并要求对随后三天的空气质量测量值进行预测。\n",
    "\n",
    "具体来说，每小时为多个站点提供为期八天的天气观测信息，例如温度，压力，风速和风向。目的是预测未来3天在多个地点的空气质量测量结果。预测的交货时间不连续；相反，必须在72小时的预测期内预测特定的交货时间。他们是：\n",
    "+1, +2, +3, +4, +5, +10, +17, +24, +48, +72\n",
    "此外，该数据集被分为不相交但连续的数据块，其中八天的数据紧随其后的三天需要预测。\n",
    "\n",
    "并非所有观察值在所有站点或块均可用，并且并非所有输出变量在所有站点或块均可用。在丢失的数据中，有很大一部分必须加以解决。\n",
    "\n",
    "该数据集被用作2012年Kaggle网站上的短期机器学习竞赛（或hackathon）的基础。\n",
    "\n",
    "比赛提交的内容是根据参与者保留的真实观察结果进行评估的，并使用平均绝对误差（MAE）进行评分。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 模型评估\n",
    "在评估朴素的预测方法之前，我们必须开发一种测试工具。\n",
    "\n",
    "这至少包括如何准备数据以及如何评估预测。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 加载数据集\n",
    "第一步是下载数据集并将其加载到内存中。\n",
    "\n",
    "我们的重点是包含训练数据集的“ TrainingData.csv ”文件，特别是大块中的数据，其中每个大块是连续八天的观测值和目标变量。\n",
    "\n",
    "我们可以使用Pandas read_csv（）函数将数据文件加载到内存中，并在第0行指定标题行。\n",
    "\n",
    "我们可以按“ chunkID”变量（列索引1）对数据进行分组。\n",
    "\n",
    "下面定义了一个名为to_chunks（）的函数，该函数接受已加载数据的NumPy数组，并将chunk_id的字典返回给该块的行。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Total Chunks: 208\n"
     ]
    }
   ],
   "source": [
    "# load data and split into chunks\n",
    "from numpy import unique\n",
    "from pandas import read_csv\n",
    " \n",
    "# split the dataset by 'chunkID', return a dict of id to rows\n",
    "def to_chunks(values, chunk_ix=1):\n",
    "\tchunks = dict()\n",
    "\t# get the unique chunk ids\n",
    "\tchunk_ids = unique(values[:, chunk_ix])\n",
    "\t# group rows by chunk id\n",
    "\tfor chunk_id in chunk_ids:\n",
    "\t\tselection = values[:, chunk_ix] == chunk_id\n",
    "\t\tchunks[chunk_id] = values[selection, :]\n",
    "\treturn chunks\n",
    " \n",
    "# load dataset\n",
    "dataset = read_csv('../AirQualityPrediction/TrainingData.csv', header=0)\n",
    "# group data by chunks\n",
    "values = dataset.values\n",
    "chunks = to_chunks(values)\n",
    "print('Total Chunks: %d' % len(chunks))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 资料准备\n",
    "现在我们知道了如何加载数据并将其拆分为大块，我们可以将其分为训练数据集和测试数据集。\n",
    "\n",
    "尽管每个块内的实际观察次数可能相差很大，但是每个块的时间间隔为每小时观察八天。\n",
    "\n",
    "我们可以将每个块分为观察的前五天进行训练，最后三天进行测试。\n",
    "\n",
    "每个观察值都有一个称为“ position_within_chunk ” 的行，范围从1到192（8天* 24小时）不等。因此，我们可以将此列中具有小于或等于120（5 * 24）的值的所有行作为训练数据，并将大于120的任何值作为测试数据。\n",
    "\n",
    "此外，在训练或测试拆分中没有任何观察结果的任何块都可以丢弃为不可行。\n",
    "\n",
    "当使用朴素模型时，我们只对目标变量感兴趣，而对输入的气象变量都不感兴趣。因此，我们可以删除输入数据，并使训练和测试数据仅包含每个块的39个目标变量，以及在块内和观察小时内的位置。\n",
    "\n",
    "下面的split_train_test（）函数实现了此行为；给定一个大块字典，它将把每个大块分解为训练和测试大块数据列表。\n",
    "\n",
    "我们不需要整个测试数据集；取而代之的是，我们只需要在三天内的特定提前期进行观察，尤其是提前期：+1, +2, +3, +4, +5, +10, +17, +24, +48, +72\n",
    "\n",
    "首先，我们可以将这些提前期放入函数中以便于参考：\n",
    "\n",
    "def get_lead_times():   \n",
    "\treturn [1, 2 ,3, 4, 5, 10, 17, 24, 48, 72]\n",
    "\n",
    "接下来，我们可以将测试数据集缩减为仅在首选交货时间处的数据。\n",
    "\n",
    "我们可以通过查看“ position_within_chunk ”列并使用提前期作为距训练数据集末尾的偏移量（例如120 + 1、120 +2等）来实现。\n",
    "\n",
    "如果在测试集中找到匹配的行，则将其保存，否则将生成一行NaN观测值。\n",
    "\n",
    "下面的函数to_forecasts（）实现了此功能，并返回一个NumPy数组，其中每个块的每个预测提前期都有一行。\n",
    "\n",
    "我们可以将所有这些结合在一起，并将数据集分为训练集和测试集，然后将结果保存到新文件中。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">dropping chunk=69: train=(0, 95), test=(28, 95)\n",
      "Train Rows: (23514, 42)\n",
      "Test Rows: (2070, 42)\n"
     ]
    }
   ],
   "source": [
    "# split data into train and test sets\n",
    "from numpy import unique\n",
    "from numpy import nan\n",
    "from numpy import array\n",
    "from numpy import savetxt\n",
    "from pandas import read_csv\n",
    " \n",
    "# split the dataset by 'chunkID', return a dict of id to rows\n",
    "def to_chunks(values, chunk_ix=1):\n",
    "\tchunks = dict()\n",
    "\t# get the unique chunk ids\n",
    "\tchunk_ids = unique(values[:, chunk_ix])\n",
    "\t# group rows by chunk id\n",
    "\tfor chunk_id in chunk_ids:\n",
    "\t\tselection = values[:, chunk_ix] == chunk_id\n",
    "\t\tchunks[chunk_id] = values[selection, :]\n",
    "\treturn chunks\n",
    " \n",
    "# split each chunk into train/test sets\n",
    "def split_train_test(chunks, row_in_chunk_ix=2):\n",
    "\ttrain, test = list(), list()\n",
    "\t# first 5 days of hourly observations for train\n",
    "\tcut_point = 5 * 24\n",
    "\t# enumerate chunks\n",
    "\tfor k,rows in chunks.items():\n",
    "\t\t# split chunk rows by 'position_within_chunk'\n",
    "\t\ttrain_rows = rows[rows[:,row_in_chunk_ix] <= cut_point, :]\n",
    "\t\ttest_rows = rows[rows[:,row_in_chunk_ix] > cut_point, :]\n",
    "\t\tif len(train_rows) == 0 or len(test_rows) == 0:\n",
    "\t\t\tprint('>dropping chunk=%d: train=%s, test=%s' % (k, train_rows.shape, test_rows.shape))\n",
    "\t\t\tcontinue\n",
    "\t\t# store with chunk id, position in chunk, hour and all targets\n",
    "\t\tindices = [1,2,5] + [x for x in range(56,train_rows.shape[1])]\n",
    "\t\ttrain.append(train_rows[:, indices])\n",
    "\t\ttest.append(test_rows[:, indices])\n",
    "\treturn train, test\n",
    " \n",
    "# return a list of relative forecast lead times\n",
    "def get_lead_times():\n",
    "\treturn [1, 2 ,3, 4, 5, 10, 17, 24, 48, 72]\n",
    " \n",
    "# convert the rows in a test chunk to forecasts\n",
    "def to_forecasts(test_chunks, row_in_chunk_ix=1):\n",
    "\t# get lead times\n",
    "\tlead_times = get_lead_times()\n",
    "\t# first 5 days of hourly observations for train\n",
    "\tcut_point = 5 * 24\n",
    "\tforecasts = list()\n",
    "\t# enumerate each chunk\n",
    "\tfor rows in test_chunks:\n",
    "\t\tchunk_id = rows[0, 0]\n",
    "\t\t# enumerate each lead time\n",
    "\t\tfor tau in lead_times:\n",
    "\t\t\t# determine the row in chunk we want for the lead time\n",
    "\t\t\toffset = cut_point + tau\n",
    "\t\t\t# retrieve data for the lead time using row number in chunk\n",
    "\t\t\trow_for_tau = rows[rows[:,row_in_chunk_ix]==offset, :]\n",
    "\t\t\t# check if we have data\n",
    "\t\t\tif len(row_for_tau) == 0:\n",
    "\t\t\t\t# create a mock row [chunk, position, hour] + [nan...]\n",
    "\t\t\t\trow = [chunk_id, offset, nan] + [nan for _ in range(39)]\n",
    "\t\t\t\tforecasts.append(row)\n",
    "\t\t\telse:\n",
    "\t\t\t\t# store the forecast row\n",
    "\t\t\t\tforecasts.append(row_for_tau[0])\n",
    "\treturn array(forecasts)\n",
    " \n",
    "# load dataset\n",
    "dataset = read_csv('../AirQualityPrediction/TrainingData.csv', header=0)\n",
    "# group data by chunks\n",
    "values = dataset.values\n",
    "chunks = to_chunks(values)\n",
    "# split into train/test\n",
    "train, test = split_train_test(chunks)\n",
    "# flatten training chunks to rows\n",
    "train_rows = array([row for rows in train for row in rows])\n",
    "# print(train_rows.shape)\n",
    "print('Train Rows: %s' % str(train_rows.shape))\n",
    "# reduce train to forecast lead times only\n",
    "test_rows = to_forecasts(test)\n",
    "print('Test Rows: %s' % str(test_rows.shape))\n",
    "# save datasets\n",
    "savetxt('naive_train.csv', train_rows, delimiter=',')\n",
    "savetxt('naive_test.csv', test_rows, delimiter=',')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "运行示例第一条注释，因为数据不足，从数据集中删除了块69。\n",
    "\n",
    "然后我们可以看到，每个训练集和测试集中都有42列，其中一列表示块ID，块内的位置，一天中的小时以及39个训练变量。\n",
    "\n",
    "我们还可以看到测试数据集的版本大大缩小，仅在预测提前期有行。\n",
    "\n",
    "新的训练和测试数据集分别保存在“ naive_train.csv ”和“ naive_test.csv ”文件中。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 预测评估\n",
    "做出预测后，需要对其进行评估。\n",
    "\n",
    "评估预测时使用更简单的格式很有帮助。例如，我们将使用[chunks] [variables] [time]的三维结构，其中variable是目标变量号（从0到38），而time是提前期索引从0到9。\n",
    "\n",
    "预计模型将以这种格式进行预测。\n",
    "\n",
    "我们还可以重组测试数据集以使该数据集进行比较。下面的prepare_test_forecasts（）函数实现了这一点。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert the test dataset in chunks to [chunk][variable][time] format\n",
    "def prepare_test_forecasts(test_chunks):\n",
    "\tpredictions = list()\n",
    "\t# enumerate chunks to forecast\n",
    "\tfor rows in test_chunks:\n",
    "\t\t# enumerate targets for chunk\n",
    "\t\tchunk_predictions = list()\n",
    "\t\tfor j in range(3, rows.shape[1]):\n",
    "\t\t\tyhat = rows[:, j]\n",
    "\t\t\tchunk_predictions.append(yhat)\n",
    "\t\tchunk_predictions = array(chunk_predictions)\n",
    "\t\tpredictions.append(chunk_predictions)\n",
    "\treturn array(predictions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "我们将使用平均绝对误差或MAE评估模型。这是比赛中使用的指标，考虑到目标变量的非高斯分布，这是明智的选择。\n",
    "\n",
    "如果提前期在测试集中不包含任何数据（例如NaN），则不会为该预测计算错误。如果提前期在测试集中确实有数据，但预测中没有数据，则观测值的全部幅度将被视为误差。最后，如果测试集具有观测值并做出了预测，则绝对差将记录为误差。\n",
    "\n",
    "该calculate_error（）函数实现这些规则，并返回错误对于一个给定的预测。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# calculate the error between an actual and predicted value\n",
    "def calculate_error(actual, predicted):\n",
    "\t# give the full actual value if predicted is nan\n",
    "\tif isnan(predicted):\n",
    "\t\treturn abs(actual)\n",
    "\t# calculate abs difference\n",
    "\treturn abs(actual - predicted)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "将所有块和所有提前期的错误求和，然后平均。\n",
    "\n",
    "将计算总体MAE，但我们还将为每个预测的提前期计算MAE。通常，这可以帮助选择模型，因为某些模型在不同的交货时间可能会有不同的表现。\n",
    "\n",
    "下面的evaluate_forecasts（）函数实现了此功能，以[chunk] [variable] [time]格式为提供的预测和期望值计算了MAE和提前期MAE 。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# evaluate a forecast in the format [chunk][variable][time]\n",
    "def evaluate_forecasts(predictions, testset):\n",
    "\tlead_times = get_lead_times()\n",
    "\ttotal_mae, times_mae = 0.0, [0.0 for _ in range(len(lead_times))]\n",
    "\ttotal_c, times_c = 0, [0 for _ in range(len(lead_times))]\n",
    "\t# enumerate test chunks\n",
    "\tfor i in range(len(test_chunks)):\n",
    "\t\t# convert to forecasts\n",
    "\t\tactual = testset[i]\n",
    "\t\tpredicted = predictions[i]\n",
    "\t\t# enumerate target variables\n",
    "\t\tfor j in range(predicted.shape[0]):\n",
    "\t\t\t# enumerate lead times\n",
    "\t\t\tfor k in range(len(lead_times)):\n",
    "\t\t\t\t# skip if actual in nan\n",
    "\t\t\t\tif isnan(actual[j, k]):\n",
    "\t\t\t\t\tcontinue\n",
    "\t\t\t\t# calculate error\n",
    "\t\t\t\terror = calculate_error(actual[j, k], predicted[j, k])\n",
    "\t\t\t\t# update statistics\n",
    "\t\t\t\ttotal_mae += error\n",
    "\t\t\t\ttimes_mae[k] += error\n",
    "\t\t\t\ttotal_c += 1\n",
    "\t\t\t\ttimes_c[k] += 1\n",
    "\t# normalize summed absolute errors\n",
    "\ttotal_mae /= total_c\n",
    "\ttimes_mae = [times_mae[i]/times_c[i] for i in range(len(times_mae))]\n",
    "\treturn total_mae, times_mae"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "一旦我们对模型进行了评估，就可以提出它。\n",
    "\n",
    "下面的summary_error（）函数首先打印模型性能的单行摘要，然后创建每个预测提前期的MAE图。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# summarize scores\n",
    "def summarize_error(name, total_mae, times_mae):\n",
    "\t# print summary\n",
    "\tlead_times = get_lead_times()\n",
    "\tformatted = ['+%d %.3f' % (lead_times[i], times_mae[i]) for i in range(len(lead_times))]\n",
    "\ts_scores = ', '.join(formatted)\n",
    "\tprint('%s: [%.3f MAE] %s' % (name, total_mae, s_scores))\n",
    "\t# plot summary\n",
    "\tpyplot.plot([str(x) for x in lead_times], times_mae, marker='.')\n",
    "\tpyplot.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "机器学习建模\n",
    "\n",
    "可以使用机器学习对问题进行建模。\n",
    "\n",
    "大多数机器学习模型并不直接支持随着时间的推移观察的概念。取而代之的是，必须将滞后观测值视为输入特征，以便进行预测。\n",
    "\n",
    "这是机器学习算法用于时间序列预测的好处。特别是，它们能够支持大量的输入功能。这些可能是一个或多个输入时间序列的滞后观察。\n",
    "\n",
    "与传统方法相比，用于时间序列预测的机器学习算法的其他一般优势包括：\n",
    "\n",
    "在变量之间的关系中支持噪声特征和噪声的能力。\n",
    "处理无关功能的能力。\n",
    "支持变量之间的复杂关系的能力。\n",
    "该数据集的挑战是需要进行多步预测。机器学习方法可用于进行多步预测的方法主要有两种：他们是：\n",
    "\n",
    "直接。开发了一个单独的模型来预测每个预测的提前期。\n",
    "递归的。开发了一个模型来进行一步预测，然后递归使用该模型，其中将先前的预测用作输入以预测后续的交货时间。\n",
    "当预测一个短的连续交付周期时，递归方法可能有意义，而当预测不连续的交付周期时，直接方法可能更有意义。鉴于我们有兴趣在三天内预测10个连续和不连续的交货时间，因此直接方法可能更适合于空气污染预测问题。\n",
    "\n",
    "数据集有39个目标变量，并且每个预测提前期，我们为每个目标变量开发一个模型。这意味着我们需要（39 * 10）390个机器学习模型。\n",
    "\n",
    "使用机器学习算法进行时间序列预测的关键是选择输入数据。我们可以考虑三个主要数据源，这些数据源可用作输入并映射到目标变量的每个预测提前期。他们是：\n",
    "\n",
    "单变量数据，例如来自预测目标变量的滞后观测值。\n",
    "多元数据，例如来自其他变量（天气和目标）的滞后观测值。\n",
    "元数据，例如有关预测日期或时间的数据。\n",
    "可以从所有块中提取数据，从而提供丰富的数据集，以学习从输入到目标预测提前期的映射。\n",
    "\n",
    "39个目标变量实际上由14个站点的12个变量组成。\n",
    "\n",
    "由于提供数据的方式，默认的建模方法是将每个可变位置都视为独立位置。可能可以按变量折叠数据，并且可以跨多个站点对变量使用相同的模型。\n",
    "\n",
    "有些变量被故意贴错标签（例如，不同数据使用具有相同标识符的变量）。但是，也许可以识别这些标签错误的变量并将其从多站点模型中排除。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 机器学习数据准备\n",
    "在探索该数据集的机器学习模型之前，我们必须以适合模型的方式准备数据。\n",
    "\n",
    "这需要两个数据准备步骤：\n",
    "\n",
    "#### 处理丢失的数据\n",
    "准备输入输出模式。\n",
    "现在，我们将重点关注39个目标变量，而忽略气象和元数据。\n",
    "\n",
    "处理丢失的数据\n",
    "每小时对39个目标变量进行5天或更少的观察。\n",
    "\n",
    "许多块没有全部五天的数据，并且所有块都不具有所有39个目标变量的数据。\n",
    "\n",
    "在块中没有目标变量数据的情况下，不需要预测。\n",
    "\n",
    "在那些大块确实有一些目标变量数据但并非全部五天有效的情况下，序列中会有缺口。在一天的观察中，这些间隔可能需要几个小时到几天，有时甚至更长。\n",
    "\n",
    "解决这些差距的三种候选策略如下：\n",
    "\n",
    "忽略差距。\n",
    "使用无间隙的数据。\n",
    "填补空白。\n",
    "我们可以忽略差距。这样做的问题是，将数据分为输入和输出时，数据将不是连续的。训练模型时，输入的数据将不一致，但可能意味着最后n个小时的数据，或者数据分散在过去n天内。这种不一致将使学习从输入到输出的映射变得非常嘈杂，对于模型而言，可能比需要的困难得多。\n",
    "\n",
    "我们只能使用没有差距的数据。这是一个不错的选择。风险在于我们可能没有太多或足够的数据来拟合模型。\n",
    "\n",
    "最后，我们可以填补空白。这称为数据插补，可以使用许多策略来填补空白。三种可能效果良好的方法包括：\n",
    "\n",
    "持续保留最后一个观测值（线性）。\n",
    "使用大块中一天中小时的中位数。\n",
    "在一天中的各个小时中使用大块中位数。\n",
    "在本教程中，我们将使用后一种方法并通过使用一天中各个块的中位数来填补空白。经过一点测试，这种方法似乎可以产生更多的训练样本和更好的模型性能。\n",
    "\n",
    "对于给定的变量，可能缺少由缺少的行定义的观察值。具体地，每个观测具有“ position_within_chunk ”。我们预计，在训练数据集每个块有120个观察，与“ positions_within_chunk ”从1到120包含边界。\n",
    "\n",
    "因此，我们可以为每个变量创建一个由120个NaN值组成的数组，使用“ position_within_chunk ”值标记该块中的所有观测值，剩下的任何内容都将标记为  NaN。然后，我们可以绘制每个变量并寻找差距。\n",
    "\n",
    "下面的variable_to_series（）函数将获取块的行和目标变量的给定列索引，并将为该变量返回一系列120个时间步长，所有可用数据均标有块中的值。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# layout a variable with breaks in the data for missing positions\n",
    "def variable_to_series(chunk_train, col_ix, n_steps=5*24):\n",
    "\t# lay out whole series\n",
    "\tdata = [nan for _ in range(n_steps)]\n",
    "\t# mark all available data\n",
    "\tfor i in range(len(chunk_train)):\n",
    "\t\t# get position in chunk\n",
    "\t\tposition = int(chunk_train[i, 1] - 1)\n",
    "\t\t# store data\n",
    "\t\tdata[position] = chunk_train[i, col_ix]\n",
    "\treturn data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们需要为每个块计算一天中小时的并行序列，我们可以使用该序列来为块中的每个变量估算特定于小时的数据。\n",
    "\n",
    "给定一系列一天中的部分时段，下面的interpolate_hours（）函数将填充一天中缺少的时段。它通过找到第一个标记的小时，然后向前计数，填写一天中的小时，然后向后执行相同的操作来做到这一点。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# interpolate series of hours (in place) in 24 hour time\n",
    "def interpolate_hours(hours):\n",
    "\t# find the first hour\n",
    "\tix = -1\n",
    "\tfor i in range(len(hours)):\n",
    "\t\tif not isnan(hours[i]):\n",
    "\t\t\tix = i\n",
    "\t\t\tbreak\n",
    "\t# fill-forward\n",
    "\thour = hours[ix]\n",
    "\tfor i in range(ix+1, len(hours)):\n",
    "\t\t# increment hour\n",
    "\t\thour += 1\n",
    "\t\t# check for a fill\n",
    "\t\tif isnan(hours[i]):\n",
    "\t\t\thours[i] = hour % 24\n",
    "\t# fill-backward\n",
    "\thour = hours[ix]\n",
    "\tfor i in range(ix-1, -1, -1):\n",
    "\t\t# decrement hour\n",
    "\t\thour -= 1\n",
    "\t\t# check for a fill\n",
    "\t\tif isnan(hours[i]):\n",
    "\t\t\thours[i] = hour % 24"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们可以调用相同的variable_to_series（）函数（上面）来创建缺少值的小时系列（列索引2），然后调用interpolate_hours（）来填补空白。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "然后，我们可以将小时数传递给任何可能使用它的归因函数。\n",
    "\n",
    "现在，我们可以尝试使用同一小时内相同系列中的值来填充大块中的缺失值。具体来说，我们将在序列中找到所有具有相同小时数的行，并计算中值。\n",
    "\n",
    "下面的impute_missing（）会接受一个块中的所有行，为该块中一天中的小时准备好的序列，以及一个缺少变量值和变量的列索引的序列。\n",
    "\n",
    "它首先检查该序列是否全部缺少数据，如果是这种情况，则由于无法进行估算，因此会立即返回。然后，它枚举序列的时间步长，当它检测到没有数据的时间步长时，它将收集同一小时内数据的序列中的所有行，并计算中位数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# impute missing data\n",
    "def impute_missing(train_chunks, rows, hours, series, col_ix):\n",
    "\t# impute missing using the median value for hour in all series\n",
    "\timputed = list()\n",
    "\tfor i in range(len(series)):\n",
    "\t\tif isnan(series[i]):\n",
    "\t\t\t# collect all rows across all chunks for the hour\n",
    "\t\t\tall_rows = list()\n",
    "\t\t\tfor rows in train_chunks:\n",
    "\t\t\t\t[all_rows.append(row) for row in rows[rows[:,2]==hours[i]]]\n",
    "\t\t\t# calculate the central tendency for target\n",
    "\t\t\tall_rows = array(all_rows)\n",
    "\t\t\t# fill with median value\n",
    "\t\t\tvalue = nanmedian(all_rows[:, col_ix])\n",
    "\t\t\tif isnan(value):\n",
    "\t\t\t\tvalue = 0.0\n",
    "\t\t\timputed.append(value)\n",
    "\t\telse:\n",
    "\t\t\timputed.append(series[i])\n",
    "\treturn imputed\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 监督学习的数据表征方式\n",
    "我们需要将每个目标变量的序列转换为具有输入和输出的行，以便适合监督的机器学习算法。\n",
    "\n",
    "具体来说，我们有一系列，例如：\n",
    "\n",
    "[1, 2, 3, 4, 5, 6, 7, 8, 9]   \n",
    "当使用2个滞后变量预测+1的提前期时，我们将序列分为输入（X）和输出（y）模式，如下所示：\n",
    "\n",
    "X,\t\t\ty   \n",
    "1, 2,\t\t3   \n",
    "2, 3,\t\t4   \n",
    "3, 4,\t\t5   \n",
    "4, 5,\t\t6   \n",
    "5, 6,\t\t7   \n",
    "6, 7,\t\t8   \n",
    "7, 8,\t\t9   \n",
    "\n",
    "这首先要求我们选择一些滞后观测值作为输入。没有正确的答案; 相反，最好测试不同的数字并查看有效的方法。\n",
    "\n",
    "然后，我们必须针对10个预测提前期中的每一个，将系列划分为有监督的学习格式。例如，用2个滞后观测值预测+24可能看起来像：\n",
    "\n",
    "X,\t\t\ty   \n",
    "1, 2,\t\t24   \n",
    "\n",
    "然后，对39个目标变量中的每一个重复此过程。\n",
    "\n",
    "然后可以跨每个块聚合为每个目标变量的每个提前期准备的模式，以提供模型的训练数据集。\n",
    "\n",
    "我们还必须准备一个测试数据集。也就是说，为每个块的每个目标变量输入数据（X），以便我们可以将其用作输入来预测测试数据集中的提前期。如果我们选择2的滞后，则测试数据集将由每个块的每个目标变量的最后两个观察值组成。非常简单。\n",
    "\n",
    "我们可以通过定义一个函数来开始，该函数将为给定的完整（插补）系列创建输入-输出模式。\n",
    "\n",
    "下面的supervised_for_lead_time（）函数将采用一个序列，将许多滞后观察值用作输入，并使用预测提前期进行预测，然后将返回从该序列中得出的输入/输出行的列表。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# created input/output patterns from a sequence\n",
    "def supervised_for_lead_time(series, n_lag, lead_time):\n",
    "\tsamples = list()\n",
    "\t# enumerate observations and create input/output patterns\n",
    "\tfor i in range(n_lag, len(series)):\n",
    "\t\tend_ix = i + (lead_time - 1)\n",
    "\t\t# check if can create a pattern\n",
    "\t\tif end_ix >= len(series):\n",
    "\t\t\tbreak\n",
    "\t\t# retrieve input and output\n",
    "\t\tstart_ix = i - n_lag\n",
    "\t\trow = series[start_ix:i] + [series[end_ix]]\n",
    "\t\tsamples.append(row)\n",
    "\treturn samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "了解这一部分很重要。\n",
    "\n",
    "我们可以测试此功能，并在一个小的测试数据集上探索不同数量的滞后变量和预测提前期。\n",
    "\n",
    "以下是一个完整的示例，该示例生成一系列20个整数，并创建一个具有两个输入滞后的序列，并预测+6提前期。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0  1  7]\n",
      " [ 1  2  8]\n",
      " [ 2  3  9]\n",
      " [ 3  4 10]\n",
      " [ 4  5 11]\n",
      " [ 5  6 12]\n",
      " [ 6  7 13]\n",
      " [ 7  8 14]\n",
      " [ 8  9 15]\n",
      " [ 9 10 16]\n",
      " [10 11 17]\n",
      " [11 12 18]\n",
      " [12 13 19]]\n"
     ]
    }
   ],
   "source": [
    "# test supervised to input/output patterns\n",
    "from numpy import array\n",
    "\n",
    "# created input/output patterns from a sequence\n",
    "def supervised_for_lead_time(series, n_lag, lead_time):\n",
    "\tdata = list()\n",
    "\t# enumerate observations and create input/output patterns\n",
    "\tfor i in range(n_lag, len(series)):\n",
    "\t\tend_ix = i + (lead_time - 1)\n",
    "\t\t# check if can create a pattern\n",
    "\t\tif end_ix >= len(series):\n",
    "\t\t\tbreak\n",
    "\t\t# retrieve input and output\n",
    "\t\tstart_ix = i - n_lag\n",
    "\t\trow = series[start_ix:i] + [series[end_ix]]\n",
    "\t\tdata.append(row)\n",
    "\treturn array(data)\n",
    "\n",
    "# define test dataset\n",
    "data = [x for x in range(20)]\n",
    "# convert to supervised format\n",
    "result = supervised_for_lead_time(data, 2, 6)\n",
    "# display result\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "运行示例将打印出结果模式，这些模式显示出滞后观测值及其相关的预测提前期。\n",
    "\n",
    "对此示例进行实验，以使其适应这种数据转换，因为它是使用机器学习算法对时间序列进行建模的关键。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在，我们可以为给定目标变量系列的每个预测提前期调用supervised_for_lead_time（）。\n",
    "\n",
    "下面的target_to_supervised（）函数实现了这一点。首先，使用上一节中开发的函数将目标变量转换为序列并进行插补。然后针对每个目标提前期创建训练样本。还创建了目标变量的测试样本。\n",
    "\n",
    "然后，针对该目标变量返回针对每个预测提前期的训练数据和测试输入数据。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create supervised learning data for each lead time for this target\n",
    "def target_to_supervised(chunks, rows, hours, col_ix, n_lag):\n",
    "\ttrain_lead_times = list()\n",
    "\t# get series\n",
    "\tseries = variable_to_series(rows, col_ix)\n",
    "\tif not has_data(series):\n",
    "\t\treturn None, [nan for _ in range(n_lag)]\n",
    "\t# impute\n",
    "\timputed = impute_missing(chunks, rows, hours, series, col_ix)\n",
    "\t# prepare test sample for chunk-variable\n",
    "\ttest_sample = array(imputed[-n_lag:])\n",
    "\t# enumerate lead times\n",
    "\tlead_times = get_lead_times()\n",
    "\tfor lead_time in lead_times:\n",
    "\t\t# make input/output data from series\n",
    "\t\ttrain_samples = supervised_for_lead_time(imputed, n_lag, lead_time)\n",
    "\t\ttrain_lead_times.append(train_samples)\n",
    "\treturn train_lead_times, test_sample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们有碎片；我们现在需要定义函数来驱动数据准备过程。\n",
    "\n",
    "此功能建立训练和测试数据集。\n",
    "\n",
    "该方法是枚举每个目标变量，并从所有块中收集每个提前期的训练数据。同时，在对测试数据集进行预测时，我们收集所需的样本作为输入。\n",
    "\n",
    "结果是训练数据集，其维度为[var] [提前期] [sample]，其中最终维度为目标变量的预测提前期的培训样本行。该函数还返回维度为[chunk] [var] [sample]的测试数据集，其中最终维度是用于对块的目标变量进行预测的输入数据。\n",
    "\n",
    "下面的data_prep（）函数实现了此行为，并采用块格式的数据和指定数量的滞后观察值用作输入。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare training [var][lead time][sample] and test [chunk][var][sample]\n",
    "def data_prep(chunks, n_lag, n_vars=39):\n",
    "\tlead_times = get_lead_times()\n",
    "\ttrain_data = [[list() for _ in range(len(lead_times))] for _ in range(n_vars)]\n",
    "\ttest_data = [[list() for _ in range(n_vars)] for _ in range(len(chunks))]\n",
    "\t# enumerate targets for chunk\n",
    "\tfor var in range(n_vars):\n",
    "\t\t# convert target number into column number\n",
    "\t\tcol_ix = 3 + var\n",
    "\t\t# enumerate chunks to forecast\n",
    "\t\tfor c_id in range(len(chunks)):\n",
    "\t\t\trows = chunks[c_id]\n",
    "\t\t\t# prepare sequence of hours for the chunk\n",
    "\t\t\thours = variable_to_series(rows, 2)\n",
    "\t\t\t# interpolate hours\n",
    "\t\t\tinterpolate_hours(hours)\n",
    "\t\t\t# check for no data\n",
    "\t\t\tif not has_data(rows[:, col_ix]):\n",
    "\t\t\t\tcontinue\n",
    "\t\t\t# convert series into training data for each lead time\n",
    "\t\t\ttrain, test_sample = target_to_supervised(chunks, rows, hours, col_ix, n_lag)\n",
    "\t\t\t# store test sample for this var-chunk\n",
    "\t\t\ttest_data[c_id][var] = test_sample\n",
    "\t\t\tif train is not None:\n",
    "\t\t\t\t# store samples per lead time\n",
    "\t\t\t\tfor lead_time in range(len(lead_times)):\n",
    "\t\t\t\t\t# add all rows to the existing list of rows\n",
    "\t\t\t\t\ttrain_data[var][lead_time].extend(train[lead_time])\n",
    "\t\t# convert all rows for each var-lead time to a numpy array\n",
    "\t\tfor lead_time in range(len(lead_times)):\n",
    "\t\t\ttrain_data[var][lead_time] = array(train_data[var][lead_time])\n",
    "\treturn array(train_data), array(test_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 完整的例子\n",
    "我们可以将所有内容结合在一起，并使用受监督的学习格式为机器学习算法准备训练和测试数据集。\n",
    "\n",
    "在预测每个预测的提前期时，我们将使用之前的12个小时的滞后观测作为输入。\n",
    "\n",
    "然后将生成的训练和测试数据集保存为二进制NumPy数组。\n",
    "\n",
    "下面列出了完整的示例。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "((39, 10), (207, 39))\n"
     ]
    }
   ],
   "source": [
    "# prepare data\n",
    "from numpy import loadtxt\n",
    "from numpy import nan\n",
    "from numpy import isnan\n",
    "from numpy import count_nonzero\n",
    "from numpy import unique\n",
    "from numpy import array\n",
    "from numpy import nanmedian\n",
    "from numpy import save\n",
    "\n",
    "# split the dataset by 'chunkID', return a list of chunks\n",
    "def to_chunks(values, chunk_ix=0):\n",
    "\tchunks = list()\n",
    "\t# get the unique chunk ids\n",
    "\tchunk_ids = unique(values[:, chunk_ix])\n",
    "\t# group rows by chunk id\n",
    "\tfor chunk_id in chunk_ids:\n",
    "\t\tselection = values[:, chunk_ix] == chunk_id\n",
    "\t\tchunks.append(values[selection, :])\n",
    "\treturn chunks\n",
    "\n",
    "# return a list of relative forecast lead times\n",
    "def get_lead_times():\n",
    "\treturn [1, 2, 3, 4, 5, 10, 17, 24, 48, 72]\n",
    "\n",
    "# interpolate series of hours (in place) in 24 hour time\n",
    "def interpolate_hours(hours):\n",
    "\t# find the first hour\n",
    "\tix = -1\n",
    "\tfor i in range(len(hours)):\n",
    "\t\tif not isnan(hours[i]):\n",
    "\t\t\tix = i\n",
    "\t\t\tbreak\n",
    "\t# fill-forward\n",
    "\thour = hours[ix]\n",
    "\tfor i in range(ix+1, len(hours)):\n",
    "\t\t# increment hour\n",
    "\t\thour += 1\n",
    "\t\t# check for a fill\n",
    "\t\tif isnan(hours[i]):\n",
    "\t\t\thours[i] = hour % 24\n",
    "\t# fill-backward\n",
    "\thour = hours[ix]\n",
    "\tfor i in range(ix-1, -1, -1):\n",
    "\t\t# decrement hour\n",
    "\t\thour -= 1\n",
    "\t\t# check for a fill\n",
    "\t\tif isnan(hours[i]):\n",
    "\t\t\thours[i] = hour % 24\n",
    "\n",
    "# return true if the array has any non-nan values\n",
    "def has_data(data):\n",
    "\treturn count_nonzero(isnan(data)) < len(data)\n",
    "\n",
    "# impute missing data\n",
    "def impute_missing(train_chunks, rows, hours, series, col_ix):\n",
    "\t# impute missing using the median value for hour in all series\n",
    "\timputed = list()\n",
    "\tfor i in range(len(series)):\n",
    "\t\tif isnan(series[i]):\n",
    "\t\t\t# collect all rows across all chunks for the hour\n",
    "\t\t\tall_rows = list()\n",
    "\t\t\tfor rows in train_chunks:\n",
    "\t\t\t\t[all_rows.append(row) for row in rows[rows[:,2]==hours[i]]]\n",
    "\t\t\t# calculate the central tendency for target\n",
    "\t\t\tall_rows = array(all_rows)\n",
    "\t\t\t# fill with median value\n",
    "\t\t\tvalue = nanmedian(all_rows[:, col_ix])\n",
    "\t\t\tif isnan(value):\n",
    "\t\t\t\tvalue = 0.0\n",
    "\t\t\timputed.append(value)\n",
    "\t\telse:\n",
    "\t\t\timputed.append(series[i])\n",
    "\treturn imputed\n",
    "\n",
    "# layout a variable with breaks in the data for missing positions\n",
    "def variable_to_series(chunk_train, col_ix, n_steps=5*24):\n",
    "\t# lay out whole series\n",
    "\tdata = [nan for _ in range(n_steps)]\n",
    "\t# mark all available data\n",
    "\tfor i in range(len(chunk_train)):\n",
    "\t\t# get position in chunk\n",
    "\t\tposition = int(chunk_train[i, 1] - 1)\n",
    "\t\t# store data\n",
    "\t\tdata[position] = chunk_train[i, col_ix]\n",
    "\treturn data\n",
    "\n",
    "# created input/output patterns from a sequence\n",
    "def supervised_for_lead_time(series, n_lag, lead_time):\n",
    "\tsamples = list()\n",
    "\t# enumerate observations and create input/output patterns\n",
    "\tfor i in range(n_lag, len(series)):\n",
    "\t\tend_ix = i + (lead_time - 1)\n",
    "\t\t# check if can create a pattern\n",
    "\t\tif end_ix >= len(series):\n",
    "\t\t\tbreak\n",
    "\t\t# retrieve input and output\n",
    "\t\tstart_ix = i - n_lag\n",
    "\t\trow = series[start_ix:i] + [series[end_ix]]\n",
    "\t\tsamples.append(row)\n",
    "\treturn samples\n",
    "\n",
    "# create supervised learning data for each lead time for this target\n",
    "def target_to_supervised(chunks, rows, hours, col_ix, n_lag):\n",
    "\ttrain_lead_times = list()\n",
    "\t# get series\n",
    "\tseries = variable_to_series(rows, col_ix)\n",
    "\tif not has_data(series):\n",
    "\t\treturn None, [nan for _ in range(n_lag)]\n",
    "\t# impute\n",
    "\timputed = impute_missing(chunks, rows, hours, series, col_ix)\n",
    "\t# prepare test sample for chunk-variable\n",
    "\ttest_sample = array(imputed[-n_lag:])\n",
    "\t# enumerate lead times\n",
    "\tlead_times = get_lead_times()\n",
    "\tfor lead_time in lead_times:\n",
    "\t\t# make input/output data from series\n",
    "\t\ttrain_samples = supervised_for_lead_time(imputed, n_lag, lead_time)\n",
    "\t\ttrain_lead_times.append(train_samples)\n",
    "\treturn train_lead_times, test_sample\n",
    "\n",
    "# prepare training [var][lead time][sample] and test [chunk][var][sample]\n",
    "def data_prep(chunks, n_lag, n_vars=39):\n",
    "\tlead_times = get_lead_times()\n",
    "\ttrain_data = [[list() for _ in range(len(lead_times))] for _ in range(n_vars)]\n",
    "\ttest_data = [[list() for _ in range(n_vars)] for _ in range(len(chunks))]\n",
    "\t# enumerate targets for chunk\n",
    "\tfor var in range(n_vars):\n",
    "\t\t# convert target number into column number\n",
    "\t\tcol_ix = 3 + var\n",
    "\t\t# enumerate chunks to forecast\n",
    "\t\tfor c_id in range(len(chunks)):\n",
    "\t\t\trows = chunks[c_id]\n",
    "\t\t\t# prepare sequence of hours for the chunk\n",
    "\t\t\thours = variable_to_series(rows, 2)\n",
    "\t\t\t# interpolate hours\n",
    "\t\t\tinterpolate_hours(hours)\n",
    "\t\t\t# check for no data\n",
    "\t\t\tif not has_data(rows[:, col_ix]):\n",
    "\t\t\t\tcontinue\n",
    "\t\t\t# convert series into training data for each lead time\n",
    "\t\t\ttrain, test_sample = target_to_supervised(chunks, rows, hours, col_ix, n_lag)\n",
    "\t\t\t# store test sample for this var-chunk\n",
    "\t\t\ttest_data[c_id][var] = test_sample\n",
    "\t\t\tif train is not None:\n",
    "\t\t\t\t# store samples per lead time\n",
    "\t\t\t\tfor lead_time in range(len(lead_times)):\n",
    "\t\t\t\t\t# add all rows to the existing list of rows\n",
    "\t\t\t\t\ttrain_data[var][lead_time].extend(train[lead_time])\n",
    "\t\t# convert all rows for each var-lead time to a numpy array\n",
    "\t\tfor lead_time in range(len(lead_times)):\n",
    "\t\t\ttrain_data[var][lead_time] = array(train_data[var][lead_time])\n",
    "\treturn array(train_data), array(test_data)\n",
    "\n",
    "# load dataset\n",
    "train = loadtxt('naive_train.csv', delimiter=',')\n",
    "test = loadtxt('naive_test.csv', delimiter=',')\n",
    "# group data by chunks\n",
    "train_chunks = to_chunks(train)\n",
    "test_chunks = to_chunks(test)\n",
    "# convert training data into supervised learning data\n",
    "n_lag = 12\n",
    "train_data, test_data = data_prep(train_chunks, n_lag)\n",
    "print(train_data.shape, test_data.shape)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save train and test sets to file\n",
    "save('supervised_train.npy', train_data)\n",
    "save('supervised_test.npy', test_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "结果是两个二进制文件，其中包含训练和测试数据集，我们可以在以下部分中加载它们，以训练和评估有关该问题的机器学习算法。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 模型评估测试工具\n",
    "在开始评估算法之前，我们需要测试工具的更多元素。\n",
    "\n",
    "首先，我们需要能够在训练数据上拟合scikit学习模型。下面的fit_model（）函数将克隆模型配置并将其适合提供的训练数据。我们将需要为每个配置的模型配备许多（360）版本，因此将大量调用此功能。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fit a single model\n",
    "def fit_model(model, X, y):\n",
    "\t# clone the model configuration\n",
    "\tlocal_model = clone(model)\n",
    "\t# fit the model\n",
    "\tlocal_model.fit(X, y)\n",
    "\treturn local_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来，我们需要为每个变量和预测提前期组合拟合模型。\n",
    "\n",
    "我们可以通过先按变量然后按交货时间枚举训练数据集来做到这一点。然后，我们可以拟合模型并将其存储在具有相同结构的列表列表中，特别是：[var] [time] [model]。\n",
    "\n",
    "下面的fit_models（）函数实现了这一点。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fit one model for each variable and each forecast lead time [var][time][model]\n",
    "def fit_models(model, train):\n",
    "\t# prepare structure for saving models\n",
    "\tmodels = [[list() for _ in range(train.shape[1])] for _ in range(train.shape[0])]\n",
    "\t# enumerate vars\n",
    "\tfor i in range(train.shape[0]):\n",
    "\t\t# enumerate lead times\n",
    "\t\tfor j in range(train.shape[1]):\n",
    "\t\t\t# get data\n",
    "\t\t\tdata = train[i, j]\n",
    "\t\t\tX, y = data[:, :-1], data[:, -1]\n",
    "\t\t\t# fit model\n",
    "\t\t\tlocal_model = fit_model(model, X, y)\n",
    "\t\t\tmodels[i][j].append(local_model)\n",
    "\treturn models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "拟合模型是较慢的部分，可以从并行化中受益，例如使用Joblib库。这留作扩展。\n",
    "\n",
    "一旦模型拟合，就可以将其用于测试数据集的预测。\n",
    "\n",
    "准备好的测试数据集首先按块组织，然后按目标变量组织。进行预测的速度很快，并且需要首先检查是否可以进行预测（我们有输入数据），如果可以，则对目标变量使用适当的模型。然后将使用这些提前期的每个直接模型来预测变量的10个预测提前期中的每一个。\n",
    "\n",
    "下面的make_predictions（）函数实现了这一点，将模型列表列表和已加载的测试数据集作为参数，并返回结构为[chunks] [var] [time]的预测数组。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# return forecasts as [chunks][var][time]\n",
    "def make_predictions(models, test):\n",
    "\tlead_times = get_lead_times()\n",
    "\tpredictions = list()\n",
    "\t# enumerate chunks\n",
    "\tfor i in range(test.shape[0]):\n",
    "\t\t# enumerate variables\n",
    "\t\tchunk_predictions = list()\n",
    "\t\tfor j in range(test.shape[1]):\n",
    "\t\t\t# get the input pattern for this chunk and target\n",
    "\t\t\tpattern = test[i,j]\n",
    "\t\t\t# assume a nan forecast\n",
    "\t\t\tforecasts = array([nan for _ in range(len(lead_times))])\n",
    "\t\t\t# check we can make a forecast\n",
    "\t\t\tif has_data(pattern):\n",
    "\t\t\t\tpattern = pattern.reshape((1, len(pattern)))\n",
    "\t\t\t\t# forecast each lead time\n",
    "\t\t\t\tforecasts = list()\n",
    "\t\t\t\tfor k in range(len(lead_times)):\n",
    "\t\t\t\t\tyhat = models[j][k][0].predict(pattern)\n",
    "\t\t\t\t\tforecasts.append(yhat[0])\n",
    "\t\t\t\tforecasts = array(forecasts)\n",
    "\t\t\t# save forecasts fore each lead time for this variable\n",
    "\t\t\tchunk_predictions.append(forecasts)\n",
    "\t\t# save forecasts for this chunk\n",
    "\t\tchunk_predictions = array(chunk_predictions)\n",
    "\t\tpredictions.append(chunk_predictions)\n",
    "\treturn array(predictions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "我们需要评估模型的列表。\n",
    "\n",
    "我们可以定义一个通用的get_models（）函数，该函数负责定义映射到已配置的scikit-learn模型对象的模型名字典。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare a list of ml models\n",
    "def get_models(models=dict()):\n",
    "\t# ...\n",
    "\treturn models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "最后，我们需要一个函数来驱动模型评估过程。\n",
    "\n",
    "给定模型字典，列举模型，首先将模型矩阵拟合到训练数据上，对测试数据集进行预测，评估预测，然后对结果进行汇总。\n",
    "\n",
    "下面的evaluate_models（）函数实现了这一点。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# evaluate a suite of models\n",
    "def evaluate_models(models, train, test, actual):\n",
    "\tfor name, model in models.items():\n",
    "\t\t# fit models\n",
    "\t\tfits = fit_models(model, train)\n",
    "\t\t# make predictions\n",
    "\t\tpredictions = make_predictions(fits, test)\n",
    "\t\t# evaluate forecast\n",
    "\t\ttotal_mae, _ = evaluate_forecasts(predictions, actual)\n",
    "\t\t# summarize forecast\n",
    "\t\tsummarize_error(name, total_mae)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "现在，我们拥有评估机器学习模型所需的一切。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 评估线性算法\n",
    "在本节中，我们将抽查一套线性机器学习算法。\n",
    "\n",
    "线性算法是那些假设输出是输入变量的线性函数的算法。这非常类似于ARIMA等经典时间序列预测模型的假设。\n",
    "\n",
    "抽查意味着对一组模型进行评估，以便大致了解可行的方法。我们对任何优于简单自回归模型AR（2）的模型都感兴趣，该模型实现了大约0.487的MAE误差。\n",
    "\n",
    "我们将使用默认配置测试八种线性机器学习算法。特别：\n",
    "\n",
    "- 线性回归\n",
    "- 套索线性回归\n",
    "- 岭回归\n",
    "- 弹性净回归\n",
    "- 胡伯回归\n",
    "- 拉索·拉斯线性回归\n",
    "- 被动攻击性回归\n",
    "- 随机梯度下降回归\n",
    "\n",
    "我们可以在get_models（）函数中定义这些模型。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare a list of ml models\n",
    "def get_models(models=dict()):\n",
    "\t# linear models\n",
    "\tmodels['lr'] = LinearRegression()\n",
    "\tmodels['lasso'] = Lasso()\n",
    "\tmodels['ridge'] = Ridge()\n",
    "\tmodels['en'] = ElasticNet()\n",
    "\tmodels['huber'] = HuberRegressor()\n",
    "\tmodels['llars'] = LassoLars()\n",
    "\tmodels['pa'] = PassiveAggressiveRegressor(max_iter=1000, tol=1e-3)\n",
    "\tmodels['sgd'] = SGDRegressor(max_iter=1000, tol=1e-3)\n",
    "\tprint('Defined %d models' % len(models))\n",
    "\treturn models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面列出了完整的代码示例。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "((39, 10), (207, 39))\n",
      "Defined 8 models\n",
      "en: 0.595 MAE\n",
      "ridge: 0.454 MAE\n",
      "huber: 0.434 MAE\n",
      "sgd: 0.457 MAE\n",
      "pa: 0.856 MAE\n",
      "lr: 0.454 MAE\n",
      "llars: 0.631 MAE\n",
      "lasso: 0.624 MAE\n"
     ]
    }
   ],
   "source": [
    "# evaluate linear algorithms\n",
    "from numpy import load\n",
    "from numpy import loadtxt\n",
    "from numpy import nan\n",
    "from numpy import isnan\n",
    "from numpy import count_nonzero\n",
    "from numpy import unique\n",
    "from numpy import array\n",
    "from sklearn.base import clone\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.linear_model import Lasso\n",
    "from sklearn.linear_model import Ridge\n",
    "from sklearn.linear_model import ElasticNet\n",
    "from sklearn.linear_model import HuberRegressor\n",
    "from sklearn.linear_model import LassoLars\n",
    "from sklearn.linear_model import PassiveAggressiveRegressor\n",
    "from sklearn.linear_model import SGDRegressor\n",
    " \n",
    "# split the dataset by 'chunkID', return a list of chunks\n",
    "def to_chunks(values, chunk_ix=0):\n",
    "\tchunks = list()\n",
    "\t# get the unique chunk ids\n",
    "\tchunk_ids = unique(values[:, chunk_ix])\n",
    "\t# group rows by chunk id\n",
    "\tfor chunk_id in chunk_ids:\n",
    "\t\tselection = values[:, chunk_ix] == chunk_id\n",
    "\t\tchunks.append(values[selection, :])\n",
    "\treturn chunks\n",
    " \n",
    "# return true if the array has any non-nan values\n",
    "def has_data(data):\n",
    "\treturn count_nonzero(isnan(data)) < len(data)\n",
    " \n",
    "# return a list of relative forecast lead times\n",
    "def get_lead_times():\n",
    "\treturn [1, 2, 3, 4, 5, 10, 17, 24, 48, 72]\n",
    " \n",
    "# fit a single model\n",
    "def fit_model(model, X, y):\n",
    "\t# clone the model configuration\n",
    "\tlocal_model = clone(model)\n",
    "\t# fit the model\n",
    "\tlocal_model.fit(X, y)\n",
    "\treturn local_model\n",
    " \n",
    "# fit one model for each variable and each forecast lead time [var][time][model]\n",
    "def fit_models(model, train):\n",
    "\t# prepare structure for saving models\n",
    "\tmodels = [[list() for _ in range(train.shape[1])] for _ in range(train.shape[0])]\n",
    "\t# enumerate vars\n",
    "\tfor i in range(train.shape[0]):\n",
    "\t\t# enumerate lead times\n",
    "\t\tfor j in range(train.shape[1]):\n",
    "\t\t\t# get data\n",
    "\t\t\tdata = train[i, j]\n",
    "\t\t\tX, y = data[:, :-1], data[:, -1]\n",
    "\t\t\t# fit model\n",
    "\t\t\tlocal_model = fit_model(model, X, y)\n",
    "\t\t\tmodels[i][j].append(local_model)\n",
    "\treturn models\n",
    " \n",
    "# return forecasts as [chunks][var][time]\n",
    "def make_predictions(models, test):\n",
    "\tlead_times = get_lead_times()\n",
    "\tpredictions = list()\n",
    "\t# enumerate chunks\n",
    "\tfor i in range(test.shape[0]):\n",
    "\t\t# enumerate variables\n",
    "\t\tchunk_predictions = list()\n",
    "\t\tfor j in range(test.shape[1]):\n",
    "\t\t\t# get the input pattern for this chunk and target\n",
    "\t\t\tpattern = test[i,j]\n",
    "\t\t\t# assume a nan forecast\n",
    "\t\t\tforecasts = array([nan for _ in range(len(lead_times))])\n",
    "\t\t\t# check we can make a forecast\n",
    "\t\t\tif has_data(pattern):\n",
    "\t\t\t\tpattern = pattern.reshape((1, len(pattern)))\n",
    "\t\t\t\t# forecast each lead time\n",
    "\t\t\t\tforecasts = list()\n",
    "\t\t\t\tfor k in range(len(lead_times)):\n",
    "\t\t\t\t\tyhat = models[j][k][0].predict(pattern)\n",
    "\t\t\t\t\tforecasts.append(yhat[0])\n",
    "\t\t\t\tforecasts = array(forecasts)\n",
    "\t\t\t# save forecasts for each lead time for this variable\n",
    "\t\t\tchunk_predictions.append(forecasts)\n",
    "\t\t# save forecasts for this chunk\n",
    "\t\tchunk_predictions = array(chunk_predictions)\n",
    "\t\tpredictions.append(chunk_predictions)\n",
    "\treturn array(predictions)\n",
    " \n",
    "# convert the test dataset in chunks to [chunk][variable][time] format\n",
    "def prepare_test_forecasts(test_chunks):\n",
    "\tpredictions = list()\n",
    "\t# enumerate chunks to forecast\n",
    "\tfor rows in test_chunks:\n",
    "\t\t# enumerate targets for chunk\n",
    "\t\tchunk_predictions = list()\n",
    "\t\tfor j in range(3, rows.shape[1]):\n",
    "\t\t\tyhat = rows[:, j]\n",
    "\t\t\tchunk_predictions.append(yhat)\n",
    "\t\tchunk_predictions = array(chunk_predictions)\n",
    "\t\tpredictions.append(chunk_predictions)\n",
    "\treturn array(predictions)\n",
    " \n",
    "# calculate the error between an actual and predicted value\n",
    "def calculate_error(actual, predicted):\n",
    "\t# give the full actual value if predicted is nan\n",
    "\tif isnan(predicted):\n",
    "\t\treturn abs(actual)\n",
    "\t# calculate abs difference\n",
    "\treturn abs(actual - predicted)\n",
    " \n",
    "# evaluate a forecast in the format [chunk][variable][time]\n",
    "def evaluate_forecasts(predictions, testset):\n",
    "\tlead_times = get_lead_times()\n",
    "\ttotal_mae, times_mae = 0.0, [0.0 for _ in range(len(lead_times))]\n",
    "\ttotal_c, times_c = 0, [0 for _ in range(len(lead_times))]\n",
    "\t# enumerate test chunks\n",
    "\tfor i in range(len(test_chunks)):\n",
    "\t\t# convert to forecasts\n",
    "\t\tactual = testset[i]\n",
    "\t\tpredicted = predictions[i]\n",
    "\t\t# enumerate target variables\n",
    "\t\tfor j in range(predicted.shape[0]):\n",
    "\t\t\t# enumerate lead times\n",
    "\t\t\tfor k in range(len(lead_times)):\n",
    "\t\t\t\t# skip if actual in nan\n",
    "\t\t\t\tif isnan(actual[j, k]):\n",
    "\t\t\t\t\tcontinue\n",
    "\t\t\t\t# calculate error\n",
    "\t\t\t\terror = calculate_error(actual[j, k], predicted[j, k])\n",
    "\t\t\t\t# update statistics\n",
    "\t\t\t\ttotal_mae += error\n",
    "\t\t\t\ttimes_mae[k] += error\n",
    "\t\t\t\ttotal_c += 1\n",
    "\t\t\t\ttimes_c[k] += 1\n",
    "\t# normalize summed absolute errors\n",
    "\ttotal_mae /= total_c\n",
    "\ttimes_mae = [times_mae[i]/times_c[i] for i in range(len(times_mae))]\n",
    "\treturn total_mae, times_mae\n",
    " \n",
    "# summarize scores\n",
    "def summarize_error(name, total_mae):\n",
    "\tprint('%s: %.3f MAE' % (name, total_mae))\n",
    " \n",
    "# prepare a list of ml models\n",
    "def get_models(models=dict()):\n",
    "\t# linear models\n",
    "\tmodels['lr'] = LinearRegression()\n",
    "\tmodels['lasso'] = Lasso()\n",
    "\tmodels['ridge'] = Ridge()\n",
    "\tmodels['en'] = ElasticNet()\n",
    "\tmodels['huber'] = HuberRegressor()\n",
    "\tmodels['llars'] = LassoLars()\n",
    "\tmodels['pa'] = PassiveAggressiveRegressor(max_iter=1000, tol=1e-3)\n",
    "\tmodels['sgd'] = SGDRegressor(max_iter=1000, tol=1e-3)\n",
    "\tprint('Defined %d models' % len(models))\n",
    "\treturn models\n",
    " \n",
    "# evaluate a suite of models\n",
    "def evaluate_models(models, train, test, actual):\n",
    "\tfor name, model in models.items():\n",
    "\t\t# fit models\n",
    "\t\tfits = fit_models(model, train)\n",
    "\t\t# make predictions\n",
    "\t\tpredictions = make_predictions(fits, test)\n",
    "\t\t# evaluate forecast\n",
    "\t\ttotal_mae, _ = evaluate_forecasts(predictions, actual)\n",
    "\t\t# summarize forecast\n",
    "\t\tsummarize_error(name, total_mae)\n",
    " \n",
    "# load supervised datasets\n",
    "train = load('supervised_train.npy', allow_pickle=True)\n",
    "test = load('supervised_test.npy', allow_pickle=True)\n",
    "print(train.shape, test.shape)\n",
    "# load test chunks for validation\n",
    "testset = loadtxt('naive_test.csv', delimiter=',')\n",
    "test_chunks = to_chunks(testset)\n",
    "actual = prepare_test_forecasts(test_chunks)\n",
    "# prepare list of models\n",
    "models = get_models()\n",
    "# evaluate models\n",
    "evaluate_models(models, train, test, actual)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "运行示例将为每种评估算法打印MAE。\n",
    "\n",
    "我们可以看到，与简单的AR模型相比，许多算法都具有技巧，MAE低于0.487。\n",
    "\n",
    "Huber回归似乎表现最佳（使用默认配置），MAE为0.434。\n",
    "\n",
    "这很有趣，因为Huber回归或带有Huber损失的鲁棒回归是一种对训练数据集中的异常值具有鲁棒性的方法。这可能表明其他方法在进行更多数据准备（例如标准化和/或离群值删除）时可能会表现更好。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 非线性算法\n",
    "我们可以使用相同的框架来评估一组非线性和集成的机器学习算法的性能。\n",
    "\n",
    "特别：\n",
    "\n",
    "#### 非线性算法\n",
    "\n",
    "- k最近邻居\n",
    "- 分类和回归树\n",
    "- 多余的树\n",
    "- 支持向量回归\n",
    "\n",
    "#### 合奏算法\n",
    "\n",
    "- Adaboost\n",
    "- 袋装决策树\n",
    "- 随机森林\n",
    "- 多余的树木\n",
    "- 梯度提升机\n",
    "\n",
    "下面的get_models（）函数定义了这9个模型。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare a list of ml models\n",
    "def get_models(models=dict()):\n",
    "\t# non-linear models\n",
    "\tmodels['knn'] = KNeighborsRegressor(n_neighbors=7)\n",
    "\tmodels['cart'] = DecisionTreeRegressor()\n",
    "\tmodels['extra'] = ExtraTreeRegressor()\n",
    "\tmodels['svmr'] = SVR()\n",
    "\t# # ensemble models\n",
    "\tn_trees = 100\n",
    "\tmodels['ada'] = AdaBoostRegressor(n_estimators=n_trees)\n",
    "\tmodels['bag'] = BaggingRegressor(n_estimators=n_trees)\n",
    "\tmodels['rf'] = RandomForestRegressor(n_estimators=n_trees)\n",
    "\tmodels['et'] = ExtraTreesRegressor(n_estimators=n_trees)\n",
    "\tmodels['gbm'] = GradientBoostingRegressor(n_estimators=n_trees)\n",
    "\tprint('Defined %d models' % len(models))\n",
    "\treturn models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "下面提供了完整的代码清单。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "((39, 10), (207, 39))\n",
      "Defined 9 models\n",
      "knn: 0.484 MAE\n",
      "extra: 0.629 MAE\n",
      "gbm: 0.450 MAE\n",
      "cart: 0.630 MAE\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-32-cc7b6714e63e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m    185\u001b[0m \u001b[0mmodels\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mget_models\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    186\u001b[0m \u001b[0;31m# evaluate models\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 187\u001b[0;31m \u001b[0mevaluate_models\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodels\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtrain\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtest\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mactual\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<ipython-input-32-cc7b6714e63e>\u001b[0m in \u001b[0;36mevaluate_models\u001b[0;34m(models, train, test, actual)\u001b[0m\n\u001b[1;32m    166\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmodel\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mmodels\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    167\u001b[0m                 \u001b[0;31m# fit models\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 168\u001b[0;31m                 \u001b[0mfits\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfit_models\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtrain\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    169\u001b[0m                 \u001b[0;31m# make predictions\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    170\u001b[0m                 \u001b[0mpredictions\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmake_predictions\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfits\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtest\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-32-cc7b6714e63e>\u001b[0m in \u001b[0;36mfit_models\u001b[0;34m(model, train)\u001b[0m\n\u001b[1;32m     57\u001b[0m                         \u001b[0mX\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     58\u001b[0m                         \u001b[0;31m# fit model\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 59\u001b[0;31m                         \u001b[0mlocal_model\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfit_model\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mX\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     60\u001b[0m                         \u001b[0mmodels\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlocal_model\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     61\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mmodels\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-32-cc7b6714e63e>\u001b[0m in \u001b[0;36mfit_model\u001b[0;34m(model, X, y)\u001b[0m\n\u001b[1;32m     42\u001b[0m         \u001b[0mlocal_model\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mclone\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmodel\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     43\u001b[0m         \u001b[0;31m# fit the model\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 44\u001b[0;31m         \u001b[0mlocal_model\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     45\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mlocal_model\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     46\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/likaiwei/anaconda2/envs/ipykernel_py2/lib/python2.7/site-packages/sklearn/ensemble/bagging.pyc\u001b[0m in \u001b[0;36mfit\u001b[0;34m(self, X, y, sample_weight)\u001b[0m\n\u001b[1;32m    245\u001b[0m             \u001b[0mReturns\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    246\u001b[0m         \"\"\"\n\u001b[0;32m--> 247\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_fit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmax_samples\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msample_weight\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msample_weight\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    248\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    249\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_fit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mX\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmax_samples\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmax_depth\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msample_weight\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/likaiwei/anaconda2/envs/ipykernel_py2/lib/python2.7/site-packages/sklearn/ensemble/bagging.pyc\u001b[0m in \u001b[0;36m_fit\u001b[0;34m(self, X, y, max_samples, max_depth, sample_weight)\u001b[0m\n\u001b[1;32m    373\u001b[0m                 \u001b[0mtotal_n_estimators\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    374\u001b[0m                 verbose=self.verbose)\n\u001b[0;32m--> 375\u001b[0;31m             for i in range(n_jobs))\n\u001b[0m\u001b[1;32m    376\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    377\u001b[0m         \u001b[0;31m# Reduce\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/likaiwei/anaconda2/envs/ipykernel_py2/lib/python2.7/site-packages/sklearn/externals/joblib/parallel.pyc\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, iterable)\u001b[0m\n\u001b[1;32m    777\u001b[0m             \u001b[0;31m# was dispatched. In particular this covers the edge\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    778\u001b[0m             \u001b[0;31m# case of Parallel used with an exhausted iterator.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 779\u001b[0;31m             \u001b[0;32mwhile\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdispatch_one_batch\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0miterator\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    780\u001b[0m                 \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_iterating\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    781\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/likaiwei/anaconda2/envs/ipykernel_py2/lib/python2.7/site-packages/sklearn/externals/joblib/parallel.pyc\u001b[0m in \u001b[0;36mdispatch_one_batch\u001b[0;34m(self, iterator)\u001b[0m\n\u001b[1;32m    623\u001b[0m                 \u001b[0;32mreturn\u001b[0m \u001b[0mFalse\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    624\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 625\u001b[0;31m                 \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_dispatch\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtasks\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    626\u001b[0m                 \u001b[0;32mreturn\u001b[0m \u001b[0mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    627\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/likaiwei/anaconda2/envs/ipykernel_py2/lib/python2.7/site-packages/sklearn/externals/joblib/parallel.pyc\u001b[0m in \u001b[0;36m_dispatch\u001b[0;34m(self, batch)\u001b[0m\n\u001b[1;32m    586\u001b[0m         \u001b[0mdispatch_timestamp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtime\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtime\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    587\u001b[0m         \u001b[0mcb\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mBatchCompletionCallBack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdispatch_timestamp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbatch\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 588\u001b[0;31m         \u001b[0mjob\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_backend\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mapply_async\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbatch\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcallback\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    589\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_jobs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mjob\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    590\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/likaiwei/anaconda2/envs/ipykernel_py2/lib/python2.7/site-packages/sklearn/externals/joblib/_parallel_backends.pyc\u001b[0m in \u001b[0;36mapply_async\u001b[0;34m(self, func, callback)\u001b[0m\n\u001b[1;32m    109\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mapply_async\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcallback\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    110\u001b[0m         \u001b[0;34m\"\"\"Schedule a func to be run\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 111\u001b[0;31m         \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mImmediateResult\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    112\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mcallback\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    113\u001b[0m             \u001b[0mcallback\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresult\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/likaiwei/anaconda2/envs/ipykernel_py2/lib/python2.7/site-packages/sklearn/externals/joblib/_parallel_backends.pyc\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, batch)\u001b[0m\n\u001b[1;32m    330\u001b[0m         \u001b[0;31m# Don't delay the application, to avoid keeping the input\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    331\u001b[0m         \u001b[0;31m# arguments in memory\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 332\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mresults\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbatch\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    333\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    334\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/likaiwei/anaconda2/envs/ipykernel_py2/lib/python2.7/site-packages/sklearn/externals/joblib/parallel.pyc\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    129\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    130\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m__call__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 131\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkwargs\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    132\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    133\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m__len__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/likaiwei/anaconda2/envs/ipykernel_py2/lib/python2.7/site-packages/sklearn/ensemble/bagging.pyc\u001b[0m in \u001b[0;36m_parallel_build_estimators\u001b[0;34m(n_estimators, ensemble, X, y, sample_weight, seeds, total_n_estimators, verbose)\u001b[0m\n\u001b[1;32m    110\u001b[0m                 \u001b[0mcurr_sample_weight\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mnot_indices_mask\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    111\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 112\u001b[0;31m             \u001b[0mestimator\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfeatures\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msample_weight\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcurr_sample_weight\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    113\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    114\u001b[0m         \u001b[0;31m# Draw samples, using a mask, and then fit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/likaiwei/anaconda2/envs/ipykernel_py2/lib/python2.7/site-packages/sklearn/tree/tree.pyc\u001b[0m in \u001b[0;36mfit\u001b[0;34m(self, X, y, sample_weight, check_input, X_idx_sorted)\u001b[0m\n\u001b[1;32m   1122\u001b[0m             \u001b[0msample_weight\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msample_weight\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1123\u001b[0m             \u001b[0mcheck_input\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcheck_input\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1124\u001b[0;31m             X_idx_sorted=X_idx_sorted)\n\u001b[0m\u001b[1;32m   1125\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1126\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Users/likaiwei/anaconda2/envs/ipykernel_py2/lib/python2.7/site-packages/sklearn/tree/tree.pyc\u001b[0m in \u001b[0;36mfit\u001b[0;34m(self, X, y, sample_weight, check_input, X_idx_sorted)\u001b[0m\n\u001b[1;32m    360\u001b[0m                                            min_impurity_split)\n\u001b[1;32m    361\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 362\u001b[0;31m         \u001b[0mbuilder\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuild\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtree_\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mX\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msample_weight\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mX_idx_sorted\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    363\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    364\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mn_outputs_\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# spot check nonlinear algorithms\n",
    "from numpy import load\n",
    "from numpy import loadtxt\n",
    "from numpy import nan\n",
    "from numpy import isnan\n",
    "from numpy import count_nonzero\n",
    "from numpy import unique\n",
    "from numpy import array\n",
    "from sklearn.base import clone\n",
    "from sklearn.neighbors import KNeighborsRegressor\n",
    "from sklearn.tree import DecisionTreeRegressor\n",
    "from sklearn.tree import ExtraTreeRegressor\n",
    "from sklearn.svm import SVR\n",
    "from sklearn.ensemble import AdaBoostRegressor\n",
    "from sklearn.ensemble import BaggingRegressor\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.ensemble import ExtraTreesRegressor\n",
    "from sklearn.ensemble import GradientBoostingRegressor\n",
    " \n",
    "# split the dataset by 'chunkID', return a list of chunks\n",
    "def to_chunks(values, chunk_ix=0):\n",
    "\tchunks = list()\n",
    "\t# get the unique chunk ids\n",
    "\tchunk_ids = unique(values[:, chunk_ix])\n",
    "\t# group rows by chunk id\n",
    "\tfor chunk_id in chunk_ids:\n",
    "\t\tselection = values[:, chunk_ix] == chunk_id\n",
    "\t\tchunks.append(values[selection, :])\n",
    "\treturn chunks\n",
    " \n",
    "# return true if the array has any non-nan values\n",
    "def has_data(data):\n",
    "\treturn count_nonzero(isnan(data)) < len(data)\n",
    " \n",
    "# return a list of relative forecast lead times\n",
    "def get_lead_times():\n",
    "\treturn [1, 2, 3, 4, 5, 10, 17, 24, 48, 72]\n",
    " \n",
    "# fit a single model\n",
    "def fit_model(model, X, y):\n",
    "\t# clone the model configuration\n",
    "\tlocal_model = clone(model)\n",
    "\t# fit the model\n",
    "\tlocal_model.fit(X, y)\n",
    "\treturn local_model\n",
    " \n",
    "# fit one model for each variable and each forecast lead time [var][time][model]\n",
    "def fit_models(model, train):\n",
    "\t# prepare structure for saving models\n",
    "\tmodels = [[list() for _ in range(train.shape[1])] for _ in range(train.shape[0])]\n",
    "\t# enumerate vars\n",
    "\tfor i in range(train.shape[0]):\n",
    "\t\t# enumerate lead times\n",
    "\t\tfor j in range(train.shape[1]):\n",
    "\t\t\t# get data\n",
    "\t\t\tdata = train[i, j]\n",
    "\t\t\tX, y = data[:, :-1], data[:, -1]\n",
    "\t\t\t# fit model\n",
    "\t\t\tlocal_model = fit_model(model, X, y)\n",
    "\t\t\tmodels[i][j].append(local_model)\n",
    "\treturn models\n",
    " \n",
    "# return forecasts as [chunks][var][time]\n",
    "def make_predictions(models, test):\n",
    "\tlead_times = get_lead_times()\n",
    "\tpredictions = list()\n",
    "\t# enumerate chunks\n",
    "\tfor i in range(test.shape[0]):\n",
    "\t\t# enumerate variables\n",
    "\t\tchunk_predictions = list()\n",
    "\t\tfor j in range(test.shape[1]):\n",
    "\t\t\t# get the input pattern for this chunk and target\n",
    "\t\t\tpattern = test[i,j]\n",
    "\t\t\t# assume a nan forecast\n",
    "\t\t\tforecasts = array([nan for _ in range(len(lead_times))])\n",
    "\t\t\t# check we can make a forecast\n",
    "\t\t\tif has_data(pattern):\n",
    "\t\t\t\tpattern = pattern.reshape((1, len(pattern)))\n",
    "\t\t\t\t# forecast each lead time\n",
    "\t\t\t\tforecasts = list()\n",
    "\t\t\t\tfor k in range(len(lead_times)):\n",
    "\t\t\t\t\tyhat = models[j][k][0].predict(pattern)\n",
    "\t\t\t\t\tforecasts.append(yhat[0])\n",
    "\t\t\t\tforecasts = array(forecasts)\n",
    "\t\t\t# save forecasts for each lead time for this variable\n",
    "\t\t\tchunk_predictions.append(forecasts)\n",
    "\t\t# save forecasts for this chunk\n",
    "\t\tchunk_predictions = array(chunk_predictions)\n",
    "\t\tpredictions.append(chunk_predictions)\n",
    "\treturn array(predictions)\n",
    " \n",
    "# convert the test dataset in chunks to [chunk][variable][time] format\n",
    "def prepare_test_forecasts(test_chunks):\n",
    "\tpredictions = list()\n",
    "\t# enumerate chunks to forecast\n",
    "\tfor rows in test_chunks:\n",
    "\t\t# enumerate targets for chunk\n",
    "\t\tchunk_predictions = list()\n",
    "\t\tfor j in range(3, rows.shape[1]):\n",
    "\t\t\tyhat = rows[:, j]\n",
    "\t\t\tchunk_predictions.append(yhat)\n",
    "\t\tchunk_predictions = array(chunk_predictions)\n",
    "\t\tpredictions.append(chunk_predictions)\n",
    "\treturn array(predictions)\n",
    " \n",
    "# calculate the error between an actual and predicted value\n",
    "def calculate_error(actual, predicted):\n",
    "\t# give the full actual value if predicted is nan\n",
    "\tif isnan(predicted):\n",
    "\t\treturn abs(actual)\n",
    "\t# calculate abs difference\n",
    "\treturn abs(actual - predicted)\n",
    " \n",
    "# evaluate a forecast in the format [chunk][variable][time]\n",
    "def evaluate_forecasts(predictions, testset):\n",
    "\tlead_times = get_lead_times()\n",
    "\ttotal_mae, times_mae = 0.0, [0.0 for _ in range(len(lead_times))]\n",
    "\ttotal_c, times_c = 0, [0 for _ in range(len(lead_times))]\n",
    "\t# enumerate test chunks\n",
    "\tfor i in range(len(test_chunks)):\n",
    "\t\t# convert to forecasts\n",
    "\t\tactual = testset[i]\n",
    "\t\tpredicted = predictions[i]\n",
    "\t\t# enumerate target variables\n",
    "\t\tfor j in range(predicted.shape[0]):\n",
    "\t\t\t# enumerate lead times\n",
    "\t\t\tfor k in range(len(lead_times)):\n",
    "\t\t\t\t# skip if actual in nan\n",
    "\t\t\t\tif isnan(actual[j, k]):\n",
    "\t\t\t\t\tcontinue\n",
    "\t\t\t\t# calculate error\n",
    "\t\t\t\terror = calculate_error(actual[j, k], predicted[j, k])\n",
    "\t\t\t\t# update statistics\n",
    "\t\t\t\ttotal_mae += error\n",
    "\t\t\t\ttimes_mae[k] += error\n",
    "\t\t\t\ttotal_c += 1\n",
    "\t\t\t\ttimes_c[k] += 1\n",
    "\t# normalize summed absolute errors\n",
    "\ttotal_mae /= total_c\n",
    "\ttimes_mae = [times_mae[i]/times_c[i] for i in range(len(times_mae))]\n",
    "\treturn total_mae, times_mae\n",
    " \n",
    "# summarize scores\n",
    "def summarize_error(name, total_mae):\n",
    "\tprint('%s: %.3f MAE' % (name, total_mae))\n",
    " \n",
    "# prepare a list of ml models\n",
    "def get_models(models=dict()):\n",
    "\t# non-linear models\n",
    "\tmodels['knn'] = KNeighborsRegressor(n_neighbors=7)\n",
    "\tmodels['cart'] = DecisionTreeRegressor()\n",
    "\tmodels['extra'] = ExtraTreeRegressor()\n",
    "\tmodels['svmr'] = SVR()\n",
    "\t# # ensemble models\n",
    "\tn_trees = 100\n",
    "\tmodels['ada'] = AdaBoostRegressor(n_estimators=n_trees)\n",
    "\tmodels['bag'] = BaggingRegressor(n_estimators=n_trees)\n",
    "\tmodels['rf'] = RandomForestRegressor(n_estimators=n_trees)\n",
    "\tmodels['et'] = ExtraTreesRegressor(n_estimators=n_trees)\n",
    "\tmodels['gbm'] = GradientBoostingRegressor(n_estimators=n_trees)\n",
    "\tprint('Defined %d models' % len(models))\n",
    "\treturn models\n",
    " \n",
    "# evaluate a suite of models\n",
    "def evaluate_models(models, train, test, actual):\n",
    "\tfor name, model in models.items():\n",
    "\t\t# fit models\n",
    "\t\tfits = fit_models(model, train)\n",
    "\t\t# make predictions\n",
    "\t\tpredictions = make_predictions(fits, test)\n",
    "\t\t# evaluate forecast\n",
    "\t\ttotal_mae, _ = evaluate_forecasts(predictions, actual)\n",
    "\t\t# summarize forecast\n",
    "\t\tsummarize_error(name, total_mae)\n",
    " \n",
    "# load supervised datasets\n",
    "train = load('supervised_train.npy', allow_pickle=True)\n",
    "test = load('supervised_test.npy', allow_pickle=True)\n",
    "print(train.shape, test.shape)\n",
    "# load test chunks for validation\n",
    "testset = loadtxt('naive_test.csv', delimiter=',')\n",
    "test_chunks = to_chunks(testset)\n",
    "actual = prepare_test_forecasts(test_chunks)\n",
    "# prepare list of models\n",
    "models = get_models()\n",
    "# evaluate models\n",
    "evaluate_models(models, train, test, actual)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "通过运行示例，我们可以看到许多算法与自动回归算法的基线相比表现良好，尽管在上一节中没有一种算法和Huber回归一样好。\n",
    "\n",
    "支持向量回归和梯度提升机可能都值得进一步研究，以分别达到0.437和0.450的MAE。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 调滞大小\n",
    "在以前的抽查实验中，滞后观察的次数被任意固定为12。\n",
    "\n",
    "我们可以改变滞后观测值的数量，并评估对MAE的影响。某些算法可能需要更多或更少的先验观察结果，但总体趋势可能会贯穿算法。\n",
    "\n",
    "使用一系列不同数量的滞后观测值准备监督学习数据集，并在每个值上拟合和评估HuberRegressor。\n",
    "\n",
    "我尝试了以下滞后观察：\n",
    "\n",
    "[1, 3, 6, 12, 24, 36, 48]\n",
    "\n",
    "结果如下：   \n",
    "1:\t0.451   \n",
    "3:\t0.445   \n",
    "6:\t0.441   \n",
    "12:\t0.434   \n",
    "24:\t0.423   \n",
    "36:\t0.422   \n",
    "48:\t0.439   \n",
    "\n",
    "\n",
    "我们可以看到总体MAE随着滞后观测次数的增加而下降的总体趋势，至少到误差开始再次上升的时候。\n",
    "\n",
    "结果表明，至少对于HuberRegressor算法而言，36个滞后观测值可能是实现0.422的MAE的良好配置。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 扩展\n",
    "本节列出了一些扩展探索的想法。\n",
    "\n",
    "- 数据准备。探索简单的数据准备（例如标准化或统计异常值消除）是否可以改善模型性能。\n",
    "- 工程特性。探索工程特征（例如，一天中的预期小时中位数）是否可以改善模型性能\n",
    "- 气象变量。探索将滞后气象变量添加到模型中是否可以提高性能。\n",
    "- 跨站点模型。探索组合相同类型的目标变量并在站点之间重复使用模型是否可以提高性能。\n",
    "- 算法调整。探索调整一些性能更好的算法的超参数是否可以提高性能。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 摘要\n",
    "在本教程中，我们发现了如何开发机器学习模型以对空气污染数据进行多步时间序列预测。\n",
    "\n",
    "具体来说，我们了解到：\n",
    "\n",
    "- 如何估算缺失值并转换时间序列数据，以便可以通过监督学习算法对其进行建模。\n",
    "- 如何开发和评估一套用于多步时间序列预测的线性算法。\n",
    "- 如何开发和评估一套用于多步时间序列预测的非线性算法。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
